home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus 2000 #4 / Amiga Plus CD - 2000 - No. 4.iso / PowerPC / Dev / PPCRelease / Libm / vax / atan2.s < prev    next >
Encoding:
Text File  |  1993-09-23  |  6.6 KB  |  220 lines
  1. # Copyright (c) 1985 Regents of the University of California.
  2. # All rights reserved.
  3. #
  4. # Redistribution and use in source and binary forms, with or without
  5. # modification, are permitted provided that the following conditions
  6. # are met:
  7. # 1. Redistributions of source code must retain the above copyright
  8. #    notice, this list of conditions and the following disclaimer.
  9. # 2. Redistributions in binary form must reproduce the above copyright
  10. #    notice, this list of conditions and the following disclaimer in the
  11. #    documentation and/or other materials provided with the distribution.
  12. # 3. All advertising materials mentioning features or use of this software
  13. #    must display the following acknowledgement:
  14. #    This product includes software developed by the University of
  15. #    California, Berkeley and its contributors.
  16. # 4. Neither the name of the University nor the names of its contributors
  17. #    may be used to endorse or promote products derived from this software
  18. #    without specific prior written permission.
  19. #
  20. # THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  21. # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  22. # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  23. # ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  24. # FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  25. # DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  26. # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  27. # HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  28. # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  29. # OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  30. # SUCH DAMAGE.
  31. #
  32. #    @(#)atan2.s    5.4 (Berkeley) 10/9/90
  33. #
  34.     .data
  35.     .align    2
  36. _sccsid:
  37. .asciz    "@(#)atan2.s    1.2 (Berkeley) 8/21/85; 5.4 (ucb.elefunt) 10/9/90"
  38.  
  39. # ATAN2(Y,X)
  40. # RETURN ARG (X+iY)
  41. # VAX D FORMAT (56 BITS PRECISION)
  42. # CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/16/85; 
  43. #
  44. #    
  45. # Method :
  46. #    1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
  47. #    2. Reduce x to positive by (if x and y are unexceptional): 
  48. #        ARG (x+iy) = arctan(y/x)          ... if x > 0,
  49. #        ARG (x+iy) = pi - arctan[y/(-x)]   ... if x < 0,
  50. #    3. According to the integer k=4t+0.25 truncated , t=y/x, the argument 
  51. #       is further reduced to one of the following intervals and the 
  52. #       arctangent of y/x is evaluated by the corresponding formula:
  53. #
  54. #          [0,7/16]       atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
  55. #       [7/16,11/16]    atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) )
  56. #       [11/16.19/16]   atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) )
  57. #       [19/16,39/16]   atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) )
  58. #       [39/16,INF]     atan(y/x) = atan(INF) + atan( -x/y )
  59. #
  60. # Special cases:
  61. # Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y).
  62. #
  63. #    ARG( NAN , (anything) ) is NaN;
  64. #    ARG( (anything), NaN ) is NaN;
  65. #    ARG(+(anything but NaN), +-0) is +-0  ;
  66. #    ARG(-(anything but NaN), +-0) is +-PI ;
  67. #    ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2;
  68. #    ARG( +INF,+-(anything but INF and NaN) ) is +-0 ;
  69. #    ARG( -INF,+-(anything but INF and NaN) ) is +-PI;
  70. #    ARG( +INF,+-INF ) is +-PI/4 ;
  71. #    ARG( -INF,+-INF ) is +-3PI/4;
  72. #    ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2;
  73. #
  74. # Accuracy:
  75. #    atan2(y,x) returns the exact ARG(x+iy) nearly rounded. 
  76. #    
  77.     .text
  78.     .align 1
  79.     .globl    _atan2
  80. _atan2 :
  81.     .word    0x0ff4
  82.     movq    4(ap),r2        # r2 = y
  83.     movq    12(ap),r4        # r4 = x
  84.     bicw3    $0x7f,r2,r0
  85.     bicw3    $0x7f,r4,r1
  86.     cmpw    r0,$0x8000        # y is the reserved operand
  87.     jeql    resop
  88.     cmpw    r1,$0x8000        # x is the reserved operand
  89.     jeql    resop
  90.     subl2    $8,sp
  91.     bicw3    $0x7fff,r2,-4(fp)    # copy y sign bit to -4(fp)
  92.     bicw3    $0x7fff,r4,-8(fp)    # copy x sign bit to -8(fp)
  93.     cmpd    r4,$0x4080        # x = 1.0 ?    
  94.     bneq    xnot1
  95.     movq    r2,r0
  96.     bicw2    $0x8000,r0        # t = |y|
  97.     movq    r0,r2            # y = |y|
  98.     brb    begin
  99. xnot1:
  100.     bicw3    $0x807f,r2,r11        # yexp
  101.     jeql    yeq0            # if y=0 goto yeq0
  102.     bicw3    $0x807f,r4,r10        # xexp
  103.     jeql    pio2            # if x=0 goto pio2
  104.     subw2    r10,r11            # k = yexp - xexp
  105.     cmpw    r11,$0x2000        # k >= 64 (exp) ?
  106.     jgeq    pio2            # atan2 = +-pi/2
  107.     divd3    r4,r2,r0        # t = y/x  never overflow
  108.     bicw2    $0x8000,r0        # t > 0
  109.     bicw2    $0xff80,r2        # clear the exponent of y
  110.     bicw2    $0xff80,r4        # clear the exponent of x
  111.     bisw2    $0x4080,r2        # normalize y to [1,2)
  112.     bisw2    $0x4080,r4        # normalize x to [1,2)
  113.     subw2    r11,r4            # scale x so that yexp-xexp=k
  114. begin:
  115.     cmpw    r0,$0x411c        # t : 39/16
  116.     jgeq    L50
  117.     addl3    $0x180,r0,r10        # 8*t
  118.     cvtrfl    r10,r10            # [8*t] rounded to int
  119.     ashl    $-1,r10,r10        # [8*t]/2
  120.     casel    r10,$0,$4
  121. L1:    
  122.     .word    L20-L1
  123.     .word    L20-L1
  124.     .word    L30-L1
  125.     .word    L40-L1
  126.     .word    L40-L1
  127. L10:    
  128.     movq    $0xb4d9940f985e407b,r6    # Hi=.98279372324732906796d0
  129.     movq    $0x21b1879a3bc2a2fc,r8    # Lo=-.17092002525602665777d-17
  130.     subd3    r4,r2,r0        # y-x
  131.     addw2    $0x80,r0        # 2(y-x)
  132.     subd2    r4,r0            # 2(y-x)-x
  133.     addw2    $0x80,r4        # 2x    
  134.     movq    r2,r10
  135.     addw2    $0x80,r10        # 2y
  136.     addd2    r10,r2            # 3y
  137.     addd2    r4,r2            # 3y+2x
  138.     divd2    r2,r0            # (2y-3x)/(2x+3y)
  139.     brw    L60
  140. L20:    
  141.     cmpw    r0,$0x3280        # t : 2**(-28)
  142.     jlss    L80
  143.     clrq    r6            # Hi=r6=0, Lo=r8=0
  144.     clrq    r8
  145.     brw    L60
  146. L30:    
  147.     movq    $0xda7b2b0d63383fed,r6    # Hi=.46364760900080611433d0
  148.     movq    $0xf0ea17b2bf912295,r8    # Lo=.10147340032515978826d-17
  149.     movq    r2,r0
  150.     addw2    $0x80,r0        # 2y
  151.     subd2    r4,r0            # 2y-x
  152.     addw2    $0x80,r4        # 2x
  153.     addd2    r2,r4            # 2x+y
  154.     divd2    r4,r0             # (2y-x)/(2x+y)
  155.     brb    L60
  156. L50:    
  157.     movq    $0x68c2a2210fda40c9,r6    # Hi=1.5707963267948966135d1
  158.     movq    $0x06e0145c26332326,r8    # Lo=.22517417741562176079d-17
  159.     cmpw    r0,$0x5100        # y : 2**57
  160.     bgeq    L90
  161.     divd3    r2,r4,r0
  162.     bisw2    $0x8000,r0         # -x/y
  163.     brb    L60
  164. L40:    
  165.     movq    $0x68c2a2210fda4049,r6    # Hi=.78539816339744830676d0
  166.     movq    $0x06e0145c263322a6,r8    # Lo=.11258708870781088040d-17
  167.     subd3    r4,r2,r0        # y-x
  168.     addd2    r4,r2            # y+x
  169.     divd2    r2,r0            # (y-x)/(y+x)
  170. L60:    
  171.     movq    r0,r10
  172.     muld2    r0,r0
  173.     polyd    r0,$12,ptable
  174.     muld2    r10,r0
  175.     subd2    r0,r8
  176.     addd3    r8,r10,r0
  177.     addd2    r6,r0
  178. L80:    
  179.     movw    -8(fp),r2
  180.     bneq    pim
  181.     bisw2    -4(fp),r0        # return sign(y)*r0
  182.     ret
  183. L90:                    # x >= 2**25 
  184.     movq    r6,r0
  185.     brb    L80
  186. pim:
  187.     subd3    r0,$0x68c2a2210fda4149,r0    # pi-t
  188.     bisw2    -4(fp),r0
  189.     ret
  190. yeq0:
  191.     movw    -8(fp),r2        
  192.     beql    zero            # if sign(x)=1 return pi
  193.     movq    $0x68c2a2210fda4149,r0    # pi=3.1415926535897932270d1
  194.     ret
  195. zero:
  196.     clrq    r0            # return 0
  197.     ret
  198. pio2:
  199.     movq    $0x68c2a2210fda40c9,r0    # pi/2=1.5707963267948966135d1
  200.     bisw2    -4(fp),r0        # return sign(y)*pi/2
  201.     ret
  202. resop:
  203.     movq    $0x8000,r0        # propagate the reserved operand
  204.     ret
  205.     .align 2
  206. ptable:
  207.     .quad    0xb50f5ce96e7abd60
  208.     .quad    0x51e44a42c1073e02
  209.     .quad    0x3487e3289643be35
  210.     .quad    0xdb62066dffba3e54
  211.     .quad    0xcf8e2d5199abbe70
  212.     .quad    0x26f39cb884883e88
  213.     .quad    0x135117d18998be9d
  214.     .quad    0x602ce9742e883eba
  215.     .quad    0xa35ad0be8e38bee3
  216.     .quad    0xffac922249243f12
  217.     .quad    0x7f14ccccccccbf4c
  218.     .quad    0xaa8faaaaaaaa3faa
  219.     .quad    0x0000000000000000
  220.